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In deep tissue photoacoustic imaging, the spatial res¬ 
olution is inherently limited by acoustic diffraction. 
Moreover, as the ultrasound attenuation increases with 
frequency, resolution is often traded-off for penetra¬ 
tion depth. Here we report on super-resolution pho¬ 
toacoustic imaging by use of multiple speckle illumina¬ 
tion. Specifically, we show that the analysis of second- 
order fluctuations of the photoacoustic images com¬ 
bined with image deconvolution enables resolving op¬ 
tically absorbing structures beyond the acoustic diffrac¬ 
tion limit. A resolution increase of almost a factor 2 is 
demonstrated experimentally. Our method introduces 
a new framework that could potentially lead to deep 
tissue photoacoustic imaging with sub-acoustic resolu¬ 
tion. 


Light scattering prevents standard optical microscopes to ob¬ 
tain well-resolved images deep inside biological tissues. In the 
past twenty years, photoacoustic (PA) imaging has been devel¬ 
oped to overcome this limitation, by imaging optical absorption 
deep inside strongly scattering tissue with the resolution of ul¬ 
trasound [1]. PA imaging relies on the unscattered ultrasonic 
waves emitted by absorbing structures under pulsed illumina¬ 
tion via thermo-elastic stress generation. It therefore provides 
images at depth in tissue with a spatial resolution limited by 
acoustic diffraction. Ultimately, the ultrasound resolution for 
biological soft tissue is limited by the attenuation of ultrasound, 
which typically increases linearly with frequency. As a result, 
the depth-to-resolution ratio of PA imaging at depth is around 
200 in practice [1, 2]. As an illustration, axial resolution down to 
5 pm and lateral resolution down to 10 pm have been reached 
with high frequency acoustic detectors at depth up to 5 mm [3]. 

In this letter, we demonstrate that the conventional acoustic- 
diffraction limit in PA imaging may be overcome by exploiting 
PA signal fluctuations, building on the super-resolution optical 
fluctuation imaging (SOFI) technique developed for fluorescence 
microscopy [4]. SOFI is based on the idea that a higher-order sta¬ 
tistical analysis of temporal fluctuations caused by fluorescence 
blinking provides a way to resolve uncorrelated fluorophores 


within a same diffraction spot. In this work, we introduce mul¬ 
tiple optical speckle illumination as a source of fluctuations for 
super-resolution PA imaging, inspired by the principle intro¬ 
duced in optics with SOFI [4] or from derived approaches using 
speckle illumination [5]. In PA imaging, multiple speckle illu¬ 
mination was initially introduced by our group as a mean to 
palliate limited-view or highpass-filtering artefacts [6]. Here, 
we demonstrate that a second-order analysis of optical speckle- 
induced PA fluctuations also provides super-resolved PA images 
beyond the acoustic diffraction limit. 

In this work, we consider PA images reconstructed from a 
set of PA signals measured with an ultrasound array. A con¬ 
ventional backprejection algorithm is used to reconstruct the 
images, and it is assumed that the reconstructed PA quantity 
A(r) may be written as a convolution: 

A(t) = [)ij(r)xJ(r)]*(i(r) (1) 

where h is the PSF corresponding to the conventional PA imag¬ 
ing process, p a the distribution of optical absorption and I the 
optical intensity pattern (see Supplement 1, sec. l.A, for a de¬ 
tailed justification of Eq. 1). Let us now consider that the region 
of interest is successively illuminated by many different speckle 
patterns 4(r) with mean (I( r)) = Iq. The following expression 
for the mean PA image (estimated from averaging the PA images 
obtained with many realizations 4(r) of the speckle illumina¬ 
tion) 

WW=Iox[f4 [ )*K [ )1 (2) 

shows that the resolution of the reconstructed image (A) (r) is 
dictated by the spatial frequency content of h( r). Under the as¬ 
sumption that the optical speckle size is much smaller than that 
of h( r), the variance image a 2 [A] (r) for uncorrelated speckles is 
given by (see Supplement 1, sec. l.B) 

« vli*) *h 2 {r) (3) 

The variance image appears as the convolution of the squared 
object by the squared PSF, which has a higher frequency content 
than the PSF itself. As a result, the variance image is expected to 
have a higher resolution as compared to the mean image. 
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Fig. 1 . a) Experimental setup: A 5 ns laser pulse is focused on 
a rotating diffuser. The resulting speckle pattern illuminates 
a collection of absorbing beads, generating ultrasounds which 
are detected with a linear ultrasonic array, b) Typical speckle 
patterns (speckle grain size ~30 pm), c) Typical point spread 
function (PSF) of the PA imaging setup (envelope FWHM: 360 =b 
25 pm in the transverse (x) direction and 450 d= 25 pm in the axial 
(z) direction). Scale bars: 200 pm. 


Our objective was to demonstrate experimentally that the 
measurement of a two-dimensional (2-D) variance image can 
provide a super-resolution PA image of the absorption distri¬ 
bution beyond the acoustic-diffraction limit. The experimental 
setup used is shown in Fig. l.a. The beam of a nanosecond 
pulsed laser (Continuum Surelite 11-10, 532 nm wavelength, 5 ns 
pulse duration, 10 Hz repetition rate) was focused on a ground 
glass diffuser (Thorlabs, 220 grit, no significant ballistic transmis¬ 
sion). The scattered light illuminated a 2-D absorbing sample 
embedded in an agarose gel block. This phantom was located 5 
cm away from the diffuser, leading to a measured speckle grain 
size of ~ 30 pm. The absorbing sample was placed in the imag¬ 
ing plane of a linear ultrasound array (Vermon, 4 MHz center 
frequency, >60% bandwidth, 128 elements, 0.33 mm pitch), con¬ 
nected to an ultrasound scanner (Aixplorer, Supersonic Imagine, 
128-channel simultaneous acquisition at 60 MS/s). A collection 
of black polyethylene microspheres (Cospheric, 50 pm and 100 
pm in diameter) was used to fabricate phantoms with isotropic 
emitters. Estimates of the PSF h( r) were measured using isolated 
50 pm diameter microspheres, while ordered patterns to be im¬ 
aged were formed using 100 pm diameter microspheres. Fig. l.c 
shows the PSF corresponding to our imaging setup. The di¬ 
mensions of the PSF, defined as the full width at half maximum 
(FWHM) of its envelope, were 360 d= 25 pm in the transverse (x) 
direction and 450 d= 25 pm in the axial (z) direction. It is worth 
noting that the anisotropy and bipolarity of the PSF are singular 
features of limited-view and limited-bandwidth PA imaging, 
and do not have their counterpart in regular all-optical imag¬ 
ing. The measurement and analysis of the PSF h( r) is further 
described in detail in Supplement 1, sec.3. 

For each sample, a set of PA images was reconstructed for 100 
uncorrelated speckle patterns, obtained by rotating the diffuser. 



Fig. 2. a) Photograph of the absorbing sample: red and blue 
arrows indicate respectively 50 pm and 100 pm diameter beads, 
b) Mean PA image over 100 speckle realisations, mimicking 
uniform illumination, c) Square root of the variance image over 
100 speckle realisations. Scale bars: 500 pm. 


The mean and variance images were then computed on a per- 
pixel basis. As described in detail in Supplement 1, sec. 2.B, 
special care was taken to reduce sources of fluctuations other 
than the multiple speckle illumination between PA acquisitions. 
For the image reconstruction, a time-domain backprojection 
algorithm was used (see Supplement 1, sec. 2.C). Square images 
of 20 mm side were reconstructed on a grid of square pixels 
(25 pm side). To demonstrate the resolution enhancement of 
our technique, we first designed the phantom shown in Fig. 2.a. 
Pairs of 100 pm diameter beads were positioned along the z 
and x axis. The distances between beads (center to center) were: 
120 pm, 140 pm and 200 pm along the z direction (from top to 
bottom), and 250 pm, 200 pm and 160 pm along the x direction 
(from left to right). 

Fig. 2.b shows the mean PA image of the sample, obtained 
from averaging over the 100 speckle realizations, and which 
corresponds to the PA image obtained under homogenous illu¬ 
mination [6]. Fig. 2.c shows the square root of the variance of the 
same data set to make a fair comparison with the mean image 
in term of units and resolution. At this stage, no significant 
resolution enhancement can be seen from the comparison of the 
mean and variance images, for both of which the pairs of beads 
are blurred by the convolution with either h( r) (Fig. 2.b) or h 2 ( r) 
(Fig. 2.c). To demonstrate that the variance image contains sub¬ 
acoustic diffraction information thanks to the higher-frequency 
content of h 2 ( r) as compared to h( r), and to remove the peculiar 
side lobes of the PA PSF h(r) observed in Fig. 2.b and Fig. 2.c, 
image deconvolution was performed based on Eqs. 2 and 3 to re¬ 
trieve the absorption distribution p a (r). The mean and variance 
images were deconvolved respectively by h( r) and h 2 ( r). In 
an ideal noise-free situation, deconvolution of an image should 
retrieve the absorbing object with no resolution limit. However, 
in the presence of noise, spatial frequencies are accurately mea¬ 
sured up to a certain limit set by the signal-to-noise ratio (SNR). 
An inversion strategy that allow accounting for the presence 
of noise in the measurement was therefore carried out to per¬ 
form the deconvolution. Retrieving ji 2 ( r) from measurements 

of the variance image modeled as cr 2 [A] (r) = h 2 ( r) * p^(r) + £ 
(with £ accouting for the experimental noise) was carried out 
by the minimization of the following constrained least-square 
functional: 

J(x) := | \h 2 * x — o' 2 [A] 11 2 + oc\ \x\\ 2 subject to x>0 (4) 

with 11 • 11 the Euclidian norm over the image space, and oc > 0 a 
regularization parameter. The constrained minimizer pro¬ 
vides a regularized solution to the deconvolution problem, 
which in turn defines an estimation of the absorption distribu¬ 
tion ]T acc := y/x^. For comparison, the exact same approach was 
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Fig. 3. a) Mean image deconvolved by the PSF; white dashed 
lines indicates the cross-section directions, b) Square root of the 
variance image deconvolved by the squared PSF. c) Horizon¬ 
tal cross-sections, blue curve: deconvolved mean image, red 
curve: square root of deconvolved variance image, d) Vertical 
cross-sections, blue curve: deconvolved mean image, red curve: 
square root of deconvolved variance image. Scale bars: 500 pm. 
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Fig. 4. a) Mean image deconvolved by the PSF. b) Square root 
of the variance image de convolved by the squared PSF; white 
dashed lines indicates the cross-section direction, c) photograph 
of the absorbing sample, red and blue arrows indicate respec¬ 
tively 50 pm and 100 pm diameter beads, d) Cross-sections, blue 
curve: mean image, red curve: square root of variance image. 
Scale bars: 200 pm. 


applied to retrieve an estimation of ]i a from the measurement of 
the mean image. In practice, the minimization of the functional 
above required adjusting the regularization parameter oc in or¬ 
der to obtain a fair resolution v.s. noise trade-off [7, Sec. 5.6], 
and choosing Ni tera tions hi the FISTA numerical method used to 
perform the minimization (see Supplement 1, sec. l.C for further 
details on the minimization approach and the corresponding 
algorithm). 

Deconvolved images are shown in Fig. 3. As mentioned 
above, the deconvolved variance image is square rooted to re¬ 
trieve the absorption distribution The parameters of the 
deconvolution algorithm were oc = 10 5 , Aerations = 800. On 
the deconvolved mean image (Fig. 3. a), we observe that only 
absorbers separated by the largest distance in each direction are 
resolved (z direction: 200 pm, x direction: 250 pm). In compar¬ 
ison, the absorbers appear much sharper on the deconvolved 
variance image (Fig. 3.b) and almost every pair is resolved on 
this image. The only exception is the upper one, for which the 
beads are nearly touching. Cross-sections along both directions 
are shown respectively on Fig. 3.c and d. These cross-sections 
show the maximum amplitude projection of the PA image along 
a =b 100pm band around the white dashed lines plotted in Fig. 3.a. 

To further investigate the influence of the PSF anisotropy, 
a second phantom was designed, in which 100 pm diameter 
beads are equally spaced from each other with a center-to-center 
distance of about 150 pm and distributed on a quarter circle. The 
corresponding results are shown in Fig. 4. The parameters of the 
deconvolution algorithm were oc = 10 6 , Alterations = 600. Not 
a single bead can be resolved on the deconvolved mean image 
(Fig. 4.a). In contrast, the 6 beads can easily be distinguished on 
the deconvolved variance image (Fig. 4.b). The last two beads 
on the right are still distinguishable but tend to smear, according 
to the lower resolution in the transverse x direction. Because 
of the PSF anisotropy (a consequence of the limited view), we 
observe that the resolution strongly depends on the angle of the 


beads pair with respect to the probe surface. 

Speckle-induced PA fluctuation imaging is an original 
method recently proposed by our group to enhance visibility 
in PA imaging. In this study we demonstrated that fluctuation 
imaging can also be used to overcome the ultrasound diffrac¬ 
tion limit in PA tomography. We showed that this method can 
separate small absorbers unresolved under standard uniform 
illumination. This is to the best of our knowledge the first super¬ 
resolution PA imaging method which does not require optical 
focusing, and therefore is applicable beyond the ballistic regime. 
Super resolution PA microscopy was reported in previous stud¬ 
ies, but it relied on light focusing and on nonlinear absorption 
mechanisms, a situation which cannot be translated to deep 
tissue PA imaging [8, 9]. 

The common approach to obtain high resolution at depth in 
PA imaging was to date to improve the acoustic detection by em¬ 
ploying higher frequencies and increasing the detection aperture. 
Here, we demonstrated that additional efforts can be deployed 
towards the illumination, in order to induce fluctuations and 
separate discrete absorbers below the acoustic diffraction limit. 

Deconvolution was found essential to reveal the super¬ 
resolution properties with the second-order fluctuation images. 
Deconvolution has already been considered in PA imaging, as 
part of the reconstruction algorithm [10] and to compensate 
for smearing induced by the spatial impulse response of the 
finite-size detectors [11]. However, no super-resolution was yet 
demonstrated , since there were no physical mechanism extend¬ 
ing the high spatial frequency information. Our approach goes 
beyond past works by considering fluctuations in PA images as 
a signal that reveal higher spatial frequencies above the noise 
level. A gain in resolution was obtained in both transverse and 
axial directions, with a resolution enhancement of at least 1.7 in 
both directions. 

The deconvolution approach was implemented with ab- 
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sorbers of size similar to that of the absorbers used to measure 
the point spread function. As an immediate drawback, decon¬ 
volution was unable to restore the actual size of the absorbers, 
which lead on the deconvolved image to reconstructed beads 
smaller than their real size. Nonetheless, this method showed 
a very good sub-acoustic resolution performance, which was 
the objective of our work. Although the beads were expected to 
be almost equally absorbing, a difference in the reconstructed 
amplitudes was noticed. The proposed method was therefore 
not shown to be quantitative, which could be attributed to the 
non-linear deconvolution scheme. Further work should be car¬ 
ried out to investigate the possibility to retrieve quantitative 
absorption information. 

Our proof-of-concept experiments were carried out at low 
medical ultrasound frequency to simplify the controlled fabrica¬ 
tion of absorbing samples (see Supplement 1, sec. 2.A). However, 
the approach could be scaled down using high frequency de¬ 
tectors, and could also be extended to 3D PA imaging. It must 
be emphasized that speckle illumination is achievable at depth 
using a seeded laser [12], whose long coherence length still pro¬ 
vides speckle formation after several centimeters of propagation 
in scattering tissue. PA fluctuations induced by speckle illu¬ 
mination are expected to decrease as the number of speckle 
grains contained in the absorbers increases [6]. Detection of 
speckle-induced fluctuations would therefore be challenging 
deep inside tissue where the speckle grains are on the order 
of half the optical wavelength, and will require special care on 
the excitation and detection stability. In the past few years, sev¬ 
eral techniques combining shaped coherent illumination and PA 
imaging have been developed. These methods aim at focusing 
light inside tissue, possibly with sub-acoustic resolution, which 
could also lead to super-resolved PA imaging [13-15]. However, 
these methods require expensive hardware and tedious exper¬ 
imental procedures. We proposed here a very simple imaging 
technique than does not require any costly equipment and nearly 
no optical alignment. For biological applications, tissue-induced 
temporal decorrelation of speckle patterns could even be ex¬ 
ploited as a source of fluctuating illumination [16]. Alternatively, 
this super-resolution method can be extended to fluctuation of 
the absorption induced by blinking or switchable [17] contrast 
agents. 
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Super-resolution photoacoustic fluctuation imaging with multiple 
speckle illumination: supplementary material 


This document provides supplementary information to "Super-resolution photoacoustic fluctuation imaging with multiple speckle 
illumination". The first section contains theoretical informations: the conditions are given to model the photoacoustic images 
reconstruction as a convolution, the variance image is shown to result from the convolution of the squared PSF, and details on the 
deconvolution algorithm used are given. The second section provides more detailed information on the experimental methods, 
including samples preparation, measurements and processing of photoacoustic signals, and the backprojection algorithm used for the 
photoacoustic reconstruction. The last section analyze some properties of the photoacoustic point spread function. 

1. THEORY 

A. Photoacoustic imaging as a convolution process 

In its simplest form, the generation and propagation of photoacoustic (PA) pressure waves p( r, t), generated by a distribution y a ( r) of 
optical absorption illuminated by a pulsed light with intensity I(r) x f(t), may be described in a homogenous acoustic medium (with 
speed of sound c s and Griineisen coefficient T) by the following equation 

p(i,t) (SI) 

Eq.Sl shows that for a given temporal pulse shape /(f), the pressure field is linearly related to the distribution ^ fl (r) x I(r) of optical 
absorption times the optical intensity pattern. Various algorithms exist that aim at reconstructing y a (. r ) x I(r) from a set of PA 
signals measured at various points located on some boundary around the sample to be imaged [S18, S19]. In the context of our work, 
we consider a conventional back-projection algorithm based on summing values of the PA signals taken at appropriate retarded 
times [S18]. As such, the reconstructed images remain linearly related to y a ( r) x I(r). Moreover, if this linear reconstruction process 
may also be considered as translation-invariant, i.e. a spatial translation of the object simply translates the reconstructed object, then 
the reconstruction process may be written as a convolution as assumed by Eq. 1 of our work. Because the PA signals are measured on 
some boundary around the region to be imaged, the translation invariance cannot be strictly verified. However, for a field of view 
small enough such that every point that it contains is reconstructed with the same point spread function (PSF), it becomes possible to 
reasonably assume that the reconstruction process is translation-invariant. In our case, this assumption was validated by measuring 
PSF at four different locations in the field of view. The four PSF appeared identical (see sec. B below), and the results shown in the 
manuscript were independant of the particular PSF that was used to perfom the deconvolutions, thus validating our assumption that 
the reconstruction process may be written as a convolution. 

B. The variance image as a convolution with h 2 ( r) 

From A(r) = [y a (j) x f(r)] * h(r) (valid only for objects lying in a small enough field of view, see sec. A above), the variance image 
a 2 [A] (r) reads 

o 2 [A](i) = Jj C(rVOMrOMr>(r —OMr —r'0*W' 

where 

C(r',r")-(I(r , )-I(r"))-(I(r , )>-(I(r")) 

is the autocovariance of the intensity fluctuation in the speckle patterns. Under the common assumption that the speckle is wide-sense 
stationary [S20, p.67], the autocovariance function C is a function of only one variable r = r' — r". This function C(r) is sharply peaked 
around its center, and its characteristic length is by definition the speckle grain size [S20, p.72]. If the speckle grains are small enough 
compared to the PSF, C(r' — r") may be considered in the double integral above proportional to a delta function S(r r — r"), yielding 
the following convolution expression: 


^l-c 2 V 2 
dt 2 s 


a^[A]{t) oc J ^ 2 (r')/z 2 (r-r')dr' = p 2 a (r) *h 2 (t) 

This expression is strictly equivalent to that found for the second-order analysis in SOFI [S4], apart from the origin of the fluctuations. 

C. Deconvolution algorithm: FISTA 

As introduced in the main manuscript, deconvolution of the variance image was performed by minimizing the following constrained 
least-square functional: 

J(x) := | \h 2 * x — o' 2 [A] 11 2 + oc\ \x\ | 2 subject to x>0 (S2) 

From a practical viewpoint, an iterative minimization algorithm is required for the numerical evaluation of Since / is a strictly 
convex functional if oc > 0, its global minimizer is asymptotically reached by any locally convergent iteration, whatever the 
initial-guess of the algorithm. For this constrained optimization task, a natural candidate is the standard projected-gradient (PG) 
since its computational burden is very low. The PG is however rather slow to converge and the FISTA iteration [S21] that achieves 
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a faster convergence is presently considered as a popular alternative and was used in our work. In comparison to the standard 
projected-gradient algorithm [S22], the FISTA method [S21] only requires the additional storage of an auxiliary component yW. Let 
be an initial-guess and y(°) = x(°), the FISTA update for k = 0,1, • • • is: 


f x^+i) — 

| y(*+!) = x( fc+1 ) + 

with 6 > 0, P+ the projection operator over R + , and g the gradient of the penalized least-square functional above that reads 

g(x) := 2 (h 2 ^J * (h 2 * x — + ocx 


(S3) 


(S4) 


where h f stands for the adjoint of the operator h. It is worth noting that all the convolution operations that appear above can be 
computed efficiently in the Fourier domain via the fast Fourier transform (FFT) and appropriate boundary conditions [S23]. Finally, 
the global convergence of Eq. (S3) is granted provided that the constant step-size 6 is adjusted within the interval (0,0) with 

0 = 1/ max(|/z(o;)| 2 + y) 

CO 


where h(co) stands for the Fourier transform of h. 

2. EXPERIMENTAL METHODS 

A. Samples preparation 

A collection of black polyethylene microspheres (Cospheric, 50 pm and 100 pm in diameter) was used to fabricate phantoms with 
isotropic emitters. Estimates of the PSF h(r) were measured using 50 pm diameter microspheres, while ordered patterns to be imaged 
were formed using 100 pm diameter microspheres. To precisely control the relative position of the beads, melted gel was first poured 
on a mold drilled with micrometer precision (Mini Mill, Minitech). The beads were then manually placed in the molded holes of the 
solidified gel. 

B. Measurements and signal processing 

To avoid other sources of fluctuation apart from the multiple speckle illumination, appropriate care was taken to eliminate noise 
and triggering-induced temporal jitters. The intensity of the laser pulses was monitored with a photodiode to compensate for the 
laser pulse-to-pulse energy fluctuations. In addition, for each speckle illumination (i.e. each diffuser position), the PA signals were 
averaged over 25 laser pulses to improve the signal-to-noise ratio. The signals were then filtered between 1 MHz and 8 MHz (3 rd 
order butterworth filter) for noise removal outside of the array bandwidth. A set of PA images was reconstructed for 100 uncorrelated 
speckle patterns. The mean and variance images of this data set were then computed on a per-pixel basis. The same procedure was 
repeated with one single speckle realization (static diffuser). The resulting variance image was then subtracted from the previous 
variance image (with rotating diffuser). This procedure was found to correct for systematic variance noise. 

C. Backprojection algorithm 

In the context of our work, we considered a conventional back-projection algorithm based on summing values of the PA signals taken 
at appropriate retarded times. More specifically, the backprojection algorithm described in [S24, eq. (20)] was implemented while 
keeping only the first time-derivative of the processed signal, assuming point-like detectors located in the centers of the elements of 
the ultrasound array. 

3. PHOTOACOUSTIC POINT SPREAD FUNCTION 

A. Measurement 

The PSF of the imaging setup was measured by concentrating light on a single 50 pm diameter bead located in the vicinity of the 
structured 100 pm bead pattern (see Fig. 2.a). This ensured that the 50 pm microsphere was the only PA source. The diffuser was 
removed from the light path during this step. 

B. Translation-invariance of the point spread function across the field-of-view 

Following the measurement approach described above, PSFs were reconstructed for 4 different locations in the field of view to confirm 
that it can be assumed as constant, as required to model the reconstruction by a convolution process. The corresponding results are 
summarized in Fig. SI. With our current ultrasonic detection and the back-projection algorithm, we showed that this approach was 
at least valid for a 3 mm x 3 mm field-of-view. On the deconvolved images, a slight difference in the PSF shapes would results in 
additional artefacts, for instance some rebounds around the reconstructed objects that are located far from the place where the PSF is 
measured. 
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Fig. SI . a) Photograph of the absorbing sample. The red arrows indicate the locations of the 50 pm diameter beads used to measure 
the PSF of the PA system, b-e) PSF of the PA imaging setup, recorded at 4 different locations: (b) Top left, (c) Top right, (d) Bottom left, 
(e) Bottom right, f) Axial cross-sections of the 4 different PSFs (along z direction): top-left (blue), top-right (red), bottom-left (green), 
bottom-right (yellow), g) Transverse cross-sections of the 4 different PSFs (along x direction). Scale bars: 500 pm. 
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Fig. S2. a) PSF of the PA imaging system, b) Enveloppe of the PSF (modulus of the Hilbert transform), c) Squared PSF of the PA 
imaging system. Scale bars: 200 pm. 


C. Analysis 

In Fig. S2, we display the PSF (Fig. S2.a) of the PA system and its envelope (Fig. S2.b, computed using Hilbert transform), and the 
squared PSF (Fig. S2.c). The conventional resolution of the imaging system (when using uniform illumination) is given by the full 
width at half maximum (FWHM) of the enveloppe of the PSF. The FWHM measured on the data shown on Fig. S2.b was 360 db 25 pm 
in the transverse x direction and 450 =L 25 pm in the axial z direction. On Fig. S2.c, we observe that the squared PSF is sharper than the 
PSF itself, which is expected to lead to a higher measurable frequency content in the fluctuation image. The spatial frequency content 
of the PSF and the squared PSF of the PA imaging setup are shown in Fig. S3. We observe that for a given noise level, higher spatial 
frequencies are measurable with the squared PSF, hence on the variance image. 


8 


















a) 


| FT (PSF) | 


c) 


10 .log 10 1 FT(PSF) | 


'e 





b) 


E 






- 0.02 - 0.01 0 0.01 0.02 
k x (|jm -1 ) 

|FT(PSF 2 )| 



k x (|Jm -i ) 



- 0.02 - 0.01 0 0.01 0.02 
k x (|jm -1 ) 



' - 0.02 - 0.01 0 0.01 0.02 
k x (|jm -1 ) 


Fig. S3, a) Two-dimensional spatial Fourier Transform (magnitude) of the PSF. b) Two-dimensional spatial Fourier Transform 
(magnitude) of the squared PSF. c-d) Same as a-b) with logarithmic scale. 
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